Computerized tomography has been used for many years, especially in the area of medical imaging, for providing nondestructive or noninvasive techniques for generating sectional views of an object. A good description of such systems appears in U.S. Pat. No. 5,414,622, which is incorporated herein by reference. These imaging systems measure the density or metabolic activity of a section of the observed object (such as a patient's body) and produce raw acquisition data. This acquisition data consists of tomographic projections or is processed to obtain tomographic projections, which may be used in a tomographic reconstruction procedure to obtain image data associated with a sectional view of the observed object. These techniques are presently employed in medical imaging systems such as X-ray computerized tomography (CT) systems, positron emission tomography (PET) systems and Single Positron Emission Computerized Tomography (SPECT) systems. Applications in fields other than medical imaging include seismology and geophysical sciences where ground penetrating acoustic radiation is employed and reflected signals are processed to obtain tomographic projections of sub-surface structures. Other applications will be apparent to one skilled in the art.
These systems can be described mathematically as follows. A section of the object to be observed by the tomographic device is represented by an image fc (x1, x2), where x1 and x2 are the spatial dimensions of the image. In the discrete case, f is an image matrix with N1×N2 elements and each member of f includes pixel information such as, grayscale value or brightness. An estimation of f is derived by a tomographic reconstruction procedure from the tomographic projections of f, denoted Y[t,a], and defined as:Y[t,α]=Rf[t,α]+W[t,α]  (Eq. 1)
where W is an additive noise, R is the discrete Radon transform, which models the tomographic projection process, a is the angle of projection and t is the position on the line integral. The discrete Radon transform is derived from its continuous version Rc, which is equivalent to the X-ray transform in two dimensions described in S. R. Deans, The Radon Transform And Some Of Its Applications (John Wiley and Sons) (1983) and expressed as
                                          (                                          R                c                            ⁢                              f                c                                      )                    ⁢                      (                          t              ,              α                        )                          =                              ∫                          -              ∞                        ∞                    ⁢                                    ∫                              -                ∞                            ∞                        ⁢                                                            f                  c                                ⁡                                  (                                                            x                      1                                        ,                                          x                      2                                                        )                                            ⁢                              δ                ⁡                                  (                                                                                    x                        1                                            ⁢                      cos                      ⁢                                                                                          ⁢                      α                                        +                                                                  x                        2                                            ⁢                      sin                      ⁢                                                                                          ⁢                      α                                        -                    t                                    )                                            ⁢                                                          ⁢                              ⅆ                                  x                  1                                            ⁢                                                          ⁢                              ⅆ                                  x                  2                                                                                        (                  Eq          .                                          ⁢          2                )            
where fc is a square integrable function δ is the Dirac mass function where δ(x)=1 when x=0 and δ(x)=0 otherwise and a is an angle from 0 to 2π radians. FIG. 5 is a diagram showing a typical relationship between R,f t and α.
One skilled in the art will recognize that there are several ways to define the discrete Radon transform based on the continuous Radon transform. Some of those methods are described in P. Toft, The Radon Transform and Implementation, Ph.D. thesis, Department of Mathematical Modeling, Technical University of Denmark (1996). Typically, a line integral along x1 cos α+x2 sin α=t is approximated by a summation of the pixel values inside the strip t−½≦n1 cos α+n2 sin α≦t+½.
The noise W from Eq. 1 is traditionally modeled as Gaussian or Poisson noise. However, since the tomographic projections Y are often processed to incorporate various corrections, such as attenuation correction, scatter correction, resolution correction and geometric correction, the resulting noise W does not actually comply with such statistical models.
Once tomographic projections are obtained, reconstruction of the desired image data proceeds in two steps. The steps may-be performed in either order. The first of these steps is backprojection. This step entails computing the inverse radon transform R−1) of the tomographic projections Y. A discrete backprojection may be directly computed by (a) radial interpolation, using techniques well known to one skilled in the art and (b) deconvolution to amplify the high frequency components of the 1-dimensional Fourier transforms of Y in the direction of t by multiplying the 1-dimensional Fourier transforms of Y in the direction of t by a high-pass filter, such as a ramp filter.
Second, regularization if performed. It follows from Eq. 1 that computing the inverse Radon transform of the tomographic projections yields:R−1(Y)=f+R−1(W),   (Eq. 3)
Where the desired image f is contaminated by additive noise R−1(W). Because the Radon transform R is a smoothing transform, computing its inverse R−1 in the presence of additive noise W is an ill-posed inverse problem, so R−1(W) represents a large additive noise. Consequently, the regularization step is necessary to remove as much of the additive noise as possible from the final reconstructed data.
The prior art teaches two approaches to regularization, both of which suffer shortcomings. The most popular regularization technique in the prior art is known as Filtered Back-Projection (FBP). FBP and its derivatives are linear filtering techniques in the Fourier space. FBP suffers from limitations due to the fact that the vectors of the Fourier basis are not adapted to represent spatially inhomogeneous data typically found in images. The second well-known technique for regularization utilize iterative statistical models (NSM) designed to implement Expectation-Maximum (EM) and Maximum A Posteriori (MAP) estimators. These techniques are described in G. McLachlan and T. Krishnan, The EM Algorithm and Extensions (New York:Wiley) (1997) and L. A. Shepp and Y. Vardi, Maximum Likelihood Reconstruction For Emission Tomography, Transactions on Medical Imaging, 1982, at 113–122. Although ISM methods can show an improvement over FBP techniques, they suffer from several drawbacks, including long computation times, lack of theoretical understanding to characterize estimation error and convergence problems.